Solving a parabolic equation
Created: 20-04-2025. Last modified: 23-04-2025.
Let us set some global options for all code chunks in this document.
# Set seed for reproducibility
set.seed(1982)
# Set global options for all code chunks
knitr::opts_chunk$set(
# Disable messages printed by R code chunks
message = FALSE,
# Disable warnings printed by R code chunks
warning = FALSE,
# Show R code within code chunks in output
echo = TRUE,
# Include both R code and its results in output
include = TRUE,
# Evaluate R code chunks
eval = TRUE,
# Enable caching of R code chunks for faster rendering
cache = FALSE,
# Align figures in the center of the output
fig.align = "center",
# Enable retina display for high-resolution figures
retina = 2,
# Show errors in the output instead of stopping rendering
error = TRUE,
# Do not collapse code and output into a single block
collapse = FALSE
)# inla.upgrade(testing = TRUE)
# remotes::install_github("inlabru-org/inlabru", ref = "devel")
# remotes::install_github("davidbolin/rspde", ref = "devel")
# remotes::install_github("davidbolin/metricgraph", ref = "devel")
library(INLA)
library(inlabru)
library(rSPDE)
library(MetricGraph)
library(grateful)
library(plotly)We want to solve the fractional diffusion equation \[\begin{equation} \label{eq:maineq} \partial_t u+(\kappa^2-\Delta_\Gamma)^{\frac{\alpha}{2}} u=f \text { on } \Gamma \times(0, T), \quad u(0)=u_0 \text { on } \Gamma, \end{equation}\] where \(u\) satisfies the Kirchhoff vertex conditions \[\begin{equation} \label{eq:Kcond} \left\{\phi\in C(\Gamma)\;\Big|\; \forall v\in V: \sum_{e\in\mathcal{E}_v}\partial_e \phi(v)=0 \right\} \end{equation}\]
Let us start by building a graph.
graph <- metric_graph$new(perform_merges = TRUE,
tolerance = list(edge_edge = 1e-3,
vertex_vertex = 1e-3,
edge_vertex = 1e-3))
graph$plot()
Let \(\alpha\in(0,2]\) and \(U_h^\tau\) denote the sequence of approximations of the solution to the weak form of problem \(\eqref{eq:maineq}\) at each time step on a mesh indexed by \(h\). Let \(z=0\) and \(U^0_h = P_hu_0\). For \(k=0,\dots, N-1\), \(U_h^{k+1}\in V_h\) solves the following scheme \[\begin{align} \label{system:fully_discrete_scheme} \langle\delta U_h^{k+1},\phi\rangle + \mathfrak{a}(U_h^{k+1},\phi) = \langle f^{k+1},\phi\rangle ,\quad\forall\phi\in V_h, \end{align}\] where \(f^{k+1} = \displaystyle\dfrac{1}{\tau}\int_{t_k}^{t^{k+1}}f(t)dt\). The solution can be represented as \[\begin{align*} U_h^k(s) = \sum_{j=1}^{N_h}u_j^k\psi_j(s) \end{align*}\] Replacing this into \(\eqref{system:fully_discrete_scheme}\) yields the following linear system \[\begin{align*} \sum_{j=1}^{N_h}u_j^{k+1}[(\psi_j,\psi_i)_{L_2(\Gamma)}+ \tau\mathfrak{a}(\psi_j,\psi_i)] = \sum_{j=1}^{N_h}u_j^{k}(\psi_j,\psi_i)_{L_2(\Gamma)}+\tau( f^{k+1},\psi_i)_{L_2(\Gamma)} \end{align*}\] for \(i = 1,\dots, N_h\). In matrix notation, \[\begin{align} \label{diff_eq_discrete} (C+\tau A)U^{k+1} = CU^k+\tau F^{k+1}, \end{align}\] where \(C\) has entries \(C_{ij} = (\psi_j,\psi_i)_{L_2(\Gamma)}\), \(A\) has entries \(A_{ij} = \mathfrak{a}(\psi_j,\psi_i)\), \(U^k\) has entries \(u_j^k\), and \(F^k\) has entries \(( f^{k},\psi_i)_{L_2(\Gamma)}\). By multiplying both sides by \(A^{-1}\) and considering its operator-based rational approximation \(P_\ell^{-1}P_r\), we arrive at \((P_rC+\tau P_\ell)U^{k+1} = P_r(CU^k+\tau F^{k+1})\).
# Compute the FEM matrices
graph$compute_fem()
G <- graph$mesh$G
C <- graph$mesh$C
x <- graph$mesh$V[, 1]
y <- graph$mesh$V[, 2]
# Initial condition
U_0 <- 10*exp(-((x-4)^2 + (y-4)^2))
# Define the time step
time_step <- 0.1
# Define the right-hand side function
fun <- function(t) {return(sin(t)*((x-4)^2 - (y-4)^2))}
# Define the time discretization
time_seq <- seq(0,pi, by = time_step)
# Compute the right-hand side function at each time step
fun_mat <- do.call(cbind, lapply(time_seq, fun))
# Define the parameters
kappa <- 1
L <- kappa^2*C + G
alpha <- 0.8
beta <- alpha/2
op <- fractional.operators(L, beta, C, scale.factor = kappa^2, m = 1)
Pl <- op$Pl
Pr <- op$Pr
funF <- C %*% fun_mat # Precompute the LHS matrix
LHS <- Pr %*% C + time_step * Pl
# Initialize U matrix to store solution at each time step
U_mat <- matrix(NA, nrow = nrow(C), ncol = length(time_seq))
U_mat[, 1] <- U_0
# Time-stepping loop
for (k in 1:(length(time_seq) - 1)) {
RHS <- Pr %*% (C %*% U_mat[, k] + time_step * funF[, k + 1])
U_mat[, k + 1] <- as.matrix(solve(LHS, RHS))
}# Plot the initial condition
p_ini <- graph$plot_function(X = U_0,
vertex_size = 1,
type = "plotly",
edge_color = "black",
edge_width = 3,
line_color = "blue",
line_width = 3)
p_ini